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Abstract 

We consider the homogenized hnear feasibihty problem, to find an 
X on the unit sphere, satisfying n hnear inequahties ajx > 0. To solve 
this problem we consider the centers of the insphere of spherical sim- 
plices, whose facets are determined by a subset of the constraints. As 
a result we find a new combinatorial algorithm for the linear feasibility 
problem. If we allow rescaling this algorithm becomes polynomial. We 
point out that the algorithm solves as well the more general convex 
feasibility problem. Moreover numerical experiments show that the 
algorithm could be of practical interest. 

Keywords: Linear programming, convex programming, feasibility 
problem, polynomial algorithm. 
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1 Introduction 

A central problem in optimization is the linear feasibility problem F. An 
instance A of F is to find for n vectors G M"^, i = 1, . . . ,n and n real 
numbers hi, i = 1, . . . ,n 

X such that aJx > z = 1, . . . , n (1) 

or to show that such an x does not exist. 
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For an instance A of F the set Pa of feasible points 



Pa = {x eR'^ \a^x>hi, ^ = 1, . . . , n} (2) 

is a polyhedron. 

The main apphcation of F is the linear optimization problem to minimize 
a hnear function over Pa- It is well known how to solve an instance of the 
linear optimization problem by solving of one or more instances of F. Of 
the vast literature on linear optimization we mention the books 1^, 0, jTTI , 
which provide the background for this paper. 

The linear feasibility problem is a particular case of the convex feasibility 
problem. Here we want to find a point in a convex set K. Clearly the problem 
depends on the representation of K. A very general way of representing K is 
by a separation oracle |]^. For every x G M'^ this oracle either confirms that 
X E K 01 gives a hyperplane separating x from K. A survey on algorithms 
for the convex feasibility problem is 0. 

There are numerous algorithms for F. Well known algorithms are the 
phase I of the simplex algorithm, Khachian's algorithm and Karmar- 
KAR's algorithm 0, for a fairly recent survey of Karmarkar's algorithm 
see e.g. 0]. Here Khachian' s algorithm and Karmarkar's algorithm are 
distinguished by the fact that they are polynomial. The simplex algorithm is 
distinguished by the fact that it is combinatorial, i.e. it's state is completely 
determined by a subset of the constraints. 

For our purposes Khachiyan's algorithm and an algorithm given by 
Agmon, Motzkin and Schoenberg |T|, @| are of particular interest as 
they are relaxation algorithms for F which may be described by the following 
simple prototype 

Algorithm 1.1 1. Choose an arbitrary x^. 

2. Choose a "simple" convex set Ck depending on x^~^ and other data 
obtained in the course of the algorithm, such that Pa C Ck and a 
suitable point x^ E Ck- If x^ G Pa terminate- 

3- Set k := k + 1 and goto step (2-)- 

In the case of Khachian's algorithm, Ck is an ellipsoid and x^ its center. 
The next ellipsoid is constructed in the following way. Choose a constraint 
z, such that ajx^ < hi. Then the next ellipsoid C^+i is the smallest ellipsoid 
circumscribed to 

Ckn{x\ ajx > ajx^}- 
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For an instance A of F and x ^ Pa say that i E {1, ... ,n} is a 
most violated constraint, if ajx — hi = minjajx ~ hj \ j = 1, . . . ,n}. In 
the Agmon-Motzkin-Schoenberg algorithm a most violated constraint, 
(jT^k y jj,^ gg^y^ ^ /c Qj^Qggj^ ^ is t^^ie projcctlon of x'^ on the half- 

space {afx > bi}. This algorithm is a very special case of an algorithm in 
i- 

Furthermore the simplex algorithm for a problem in standard form may 
be viewed as a relaxation algorithm for the dual linear feasibility problem, 
where each iterate x'^ is determined by the requirement that it satisfies d of 
the inequalities in with equality. 

Further we remark that Khachiyan's algorithm solves the problem of 
constructing a point in a convex set given by a separation oracle. We shall 
return to this topic after the presentation of our algorithms. 

We proceed as follows: In section 2 we present a relaxation algorithm 
for the homogenized form of F and determine its properties. In section 3 
we use additional transformations of the problem to construct a polynomial 
algorithm for F. Moreover we shortly discuss its relation to the convex feasi- 
bility problem. While in the previous sections we have concentrated on the 
geometric aspects of the algorithms, we discuss in section 4 the underlying 
linear algebra and present the results of some numerical experiments. 

2 Algorithms 

In a certain way we want to combine the algorithm of Agmon-Motzkin- 
SCHOENBERG with Khachiyan's algorithm. This means that we want to 
use a circumscribed body of P4 and a point, which we may consider as a 
projection on this body. Furthermore we want this body as closely related 
to Pa as possible. A rather obvious choice is a simplex bounded by d + 1 of 
the hyperplanes {aJx = bi} and the center of its inball, as this point may 
be considered as a simultaneous projection on all half-spaces. For several 
reasons this does not work very well in Euclidean space and thus we first use 
the well known process of homogenization to transform the problem into a 
spherical problem. 

To proceed we need some more notation. We write B{c,r) = G M'^ | 
||a; — c|| < r} for the ball with center c and radius r, S{c,r) = {s G M'^ | 
||c — xjl = r} for the sphere bounding B{c, r) and 5'°'"^ = ^(0, 1) for the unit 
sphere in M"'. 

Then we define the homogenization 

$ : M'^ ^ 5^ . . . , a;rf) = (xi, . . . , X,, 1)/|| (xi, . . . , x,, 1) || . (3) 
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For X — (xi, . . . , Xd+i) e W^'^^ we denote hj x — (xi, . . . , xa) its projection 
on W^. For an instance A of F we have for the set of feasible points Pa 

$(Pa) = {xeS'^\ ajx - biXd+i >0, t=l,...,n, x^+i > 0}. 

Thus wc obtain an algorithm for the linear feasibility problem F, if we have 
an algorithm for the spherical feasibility problem F', where an instance A of 
F' is to find 

X e S"^ such that aJx > 0, i = 1, . . . , n (4) 

or to show that such an x does not exist. For an instance A of F' we denote 
again by Pa the set of its feasible points. 

To study F' it is helpful to use the notion of violation. For an instance A 
of F' and x e S''' we define the violation v{x) by 

v{x) = max{0, max{— of x | i = 1, . . . , n}}. 

Thus a point is feasible for an instance ^4 of F, if and only if its violation is 
0. 

For our algorithms it is somehow more appropriate to take a polar point 
of view, i.e. we work with the normals Oj rather than with the half spheres 
aJx > 0. Moreover we construct a sequence of points, which show that the 

instance is close to being infeasible. To be more precise we denote the origin 
of ]R''+^ by O, the convex hull, ajfine hull respectively, of a set M C W^~^^ by 
conv M, aff M respectively, and define 

Definition 2.1 Let Q — {x* e M'^+-^ \ i = 1, . . . , k} be a set of points. We 
say that Q is positively spanning ifQ is affinely independent and O e conv Q. 

lik — d-\-l then Q is a positive basis. Here we also consider the case that 
dim(affg) < 

Definition 2.2 Let Q = G W^^^ \ i = 1, . . . , k} be a set of points. We 
say that Q is nearly positively spanning if Q is affinely independent and the 
orthogonal projection O' of O on aff Q is contained in convQ- The distance 
of O' to the origin is called the deficiency of Q and denoted by def Q. 

Thus a nearly positively spanning set is positively spanning if and only if 
its deficiency is 0. Moreover for a nearly positively spanning set the deficiency 
is equal to the distance of the origin to its convex hull. 

The evident significance of positively spanning sets for F' is given by 
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Proposition 2.3 Let an instance A of F' be given and Q C {ai, . . . , a„,} be 
a positively spanning set. Then every x E Pa satisfies afx = for all E Q. 

If the cardinality of Q is + 2 then A is infeasible. In the other case Pa 
is contained in the plane {ajx = | Oi G Q} and thus the dimension of the 
instance is reduced. In this case we may change the problem slightly in a 
way that Q is transformed into a nearly positively spanning set with positive 
deficiency without changing the feasibility or infeasibility of the instance. 
This is quite analogous to cases of degeneration in Khachian's algorithm 
and will not be discussed further (but compare the remarks after Lemma 

If the points of Q are contained in S*^ then positively spanning sets can 
be nicely characterized. 

Definition 2.4 Let Q be a set of affinely independent points. The sphere 
S{C, R) is said to be touching for Q, if C E aSQ and Q C S{C, R). 

Clearly every set of affinely independent points has a unique touching 
sphere. Moreover we observe that it is easy to compute a touching sphere 
for given Q. The relation of touching spheres and nearly positively spanning 
sets is given by 

Proposition 2.5 Let Q d S'^ be affinely independent, S{C,R) be its touch- 
ing sphere. Then Q is nearly positively spanning if and only if C E convQ. 
In this case def Q = y/1 — R'^. 

For further use we remark that in case C G convQ the touching sphere 
coincides with the circumsphere of Q. By means of nearly positively spanning 
sets we may describe an algorithm for F': 

Algorithm 2.6 1. Let j G {1, . . . ,n} be arbitrary, x-^ = aj, Qi = {%}. 

2. Ifx^/\\x^\\ is feasible, then stop with feasibility. Else let m be the index 
of a most violated constraint for x^ /\\x^\\, y = x^ . 

3. If Qk U {ctm} is positively spanning, then stop the algorithm. 

4. Compute the center C of the touching sphere of Qk U {a^}- 

5. If C e conv((5fc U {am}), then let x^^^ = C, Qk+i = Qk ^ {ctm}, 
k = k + 1. Goto step 2. 
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6. Let y he the point, where the line segment yC intersects the relative 
boundary of conv [Q^ U {ctm})- Let F he a facet of conv [Q^ U {a„i}), 
which contains y, and aj be the vertex of Q^, which is not contained in 
F. Let Qk = Qk \ Wj} o,nd goto step 4- 

This algorithm shares a remarkable property with the simplex algorithm. 

Definition 2.7 An algorithm for F, F' is said to be combinatorial, if the 
iterates depend only on subsets of the set of constraints. 

Thus for an instance of a combinatorial algorithm there are only finitely 



many possible iterates. As the iterates of Algorithm |2.6| are centers of the 
touching sphere of Qk, this algorithm is combinatorial like the simplex al- 
gorithm and different from Khachiyan's or Karmarkar's algorithm. To 
prove the convergence of a combinatorial algorithm it is clearly sufficient to 
show, that no iterate can be repeated. 



Lemma 2.8 Algorithm 2.6 solves F' for feasible as well as infeasible prob- 



lems in finitely many steps. 

Proof. By construction def Qk is strictly monotonely decreasing. □ 

It is certainly desirable to have some quantitative information on the 
progress of the algorithm. This can be obtained by a small change in Algo- 



rithm 2.6: 



Algorithm 2.9 Construct the pointy in step 2 of the algorithm by choosing 
y as the point on the line x^am, which is closest to the origin. 



Lemma 2.10 Algorithm \2.^ solves the feasibility problem F' and it holds 

Proof. The y constructed in step 2 has the distance from the origin given by 
the right side of (^) and by construction the distance of x'^^^ from the origin 
is certainly not greater than the distance of y from the origin. □ 

While the estimate (^) is not linear, it is independent of the dimension 
and in particular of the size of the problem. The formula may be used to 
estimate the number of steps, that the algorithm needs to reach a given level. 
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Lemma 2.11 In Algorithm \2.f\ it holds def Qk ^ after at most t^ steps. 



Proof. It is somewhat easier to consider instead of defQ^ the quantities 
Hk = (def Qfc)"^- We have from Lemma |2.1CI| and using v{x^) > 



Vk+i > Vk 



1. 



Together with yi = 1 this gives inductively ?/|^_^ > A; + 1. 



(6) 
□ 



It is not known (and was not investigated), whether the sequence of it- 
erates in Algorithm and Algorithm [^.9| can differ for the same problem. 

r9| takes large steps, if the deficiency 



According to Lemma |2.10| Algorithm 
or the violation are large, in particular this is the case at the start of the 
algorithm. 

We don't have information on the global behavior of the algorithm beyond 



Lemma |2.11| . We may obtain some more insight in its working by analyzing 

Using an appropriate system of coordinates we 



} for 



one step of Algorithm plj . 

may assume that Qk spans the affine plane H = {(a;i,0,e) | xi 
some fixed e = def > and appropriate /. Then we have x^' = (0,0,e). 
Further we have that the constraint which is chosen in step 2 has the form 
dm = {xi,X2,e + rf) with xi G M', X2 G W^~^ and rj < — e. We find for the 
center C of the touching sphere of Qk U {dm} in step 4 



C=(0,0,e) + 



-rjt 



F2lr + v 



From this we find ||C|| ~ if l^^l -C ||a;2||. In this case the algorithm 

can only increase the cardinality of Q^, but not significantly decrease def Qk- 
Clearly the crucial point of the algorithm is whether we may expect C G 
conv {Qk U {ttm})- To see what happens, let us denote by the orthogonal 
projection of onto H . By construction we have that x^ is the orthogonal 
projection of C onto H. Finally let y G relbd(conv Qfe) be the point that 
x'' G ya'm- Now an examination of the triangle y, a'^, shows that C ^ 



If we define for 



conv {Qk U {flm}) if \\y — X \\ <C ||a^ — x \\. Clearly this may happen and in 
a single step the algorithm may make very little progress. Fortunately this 
is not the complete story. 

In the algorithm we gather information about the Oj 
xeS'^ 

M{x) = {aeS'^\ 

then we have for every iterate x^ that Oj 
all ttj satisfy for all k 



a^x > —v{x)}, 



G M{x^ 



,n. 



Consequently 



tti G 



Pi m{x^ 
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Unfortunately the sets M{x^) are not convex and this information is hard to 
use. 

If we want to repeat the bad step from above, the set M{x^) tells us that 
either for the new e', rj' holds |e' + ry'l ^ |e + ^| the angle between and 
ttmi must be large. Thus we can hope that bad steps cannot be repeated too 
often. A precise analysis of the phenomena appears to be difficult. 

3 A polynomial algorithm 

In view of the last remark in the previous section it is near at hand to use 
periodically a rescaling to increase the deficiency in order to speed up the 
convergence. While it is simple to give a transformation of the sphere which 
takes polyhedra to polyhedra and increases the deficiency of a set it ap- 
pears difficult to work out the consequences of iterating such transformations. 

Fortunately it turns out that for def {Qk) sufficiently small it is possible 
to find a transformation which increases the volume of the sets Pa- To work 
out our ideas it is necessary to recall some concepts of spherical convexity. A 
set M C 5''^ is called convex, if there exists a half-sphere which contains M 
and for every pair x,y E M there exists a great circle C such that C n M is 
connected. The intersection of convex sets is itself convex. By this we may 
define for every M G S'^ which is contained in a half-sphere the spherical 
convex hull co M as the intersection of all convex subsets of S"^ which contain 
M. Here we have chosen the notation coM to distinguish the spherical 
convex hull from the Euclidean convex hull used in the previous section. 
As every half-sphere is convex we have that the intersection of half-spheres 
is convex. We say that the intersection of finitely many half-spheres is a 
(spherical) polyhedron. An example of a spherical polyhedron is the set Pa 
of feasible points of an instance A of problem F'. 

We denote for two points x,y E the spherical distance by Z{x, y), i.e. 
cos(Z(x, y)) — x'^y. Next for x e and < p < tt 

K,{x) = {yeS''\Z{x,y)<p} 

denotes the spherical cap with center x and radius p. For closed M G S'^ 
we may now define the spherical counterparts of inradius, circumradius and 
diameter. The spherical inradius r'{M), the spherical incenter c(M) re- 
spectively, are the radius, center respectively, of a largest spherical cap 
Kr>{M){c{M)) C M. We remark that in contrast to the Euclidean case 
for convex M the spherical incenter is always well defined. The spherical 
circumradius R'{M), spherical circumcenter C(M) respectively, are the ra- 
dius, center respectively, of the smallest spherical cap Kri(^m){C{M)) Z) M. 
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Here the notions R , r' were chosen to distinguish spherical circumradius and 
inradius from their Euchdean counterparts. For M C 5'°' we denote the 
boundary of K'r'(m)(C'(^)) as circumsphere of M. For any M contained in 
a half-sphere the spherical diameter D{M) is given by 



For M C S'^, ai, 02 € S'^ satisfying afx > for all x G M, afx = for at 
least one x G M and ai ^ 02 we denote by vr — Z(ai, 02) a width of M. The 
breadth d{M) is then defined as the minimum over all widths of M. 

For M G S'^ the polar set M* is defined by 



In our context it is of importance to observe that the polar set of a 
spherical polyhedron is again a polyhedron and that inclusions are reversed, 
i.e. for M C M' we have M* D M'* . For M = {x\ . . . , x^}, x\ . . . , x'' 
linearly independent, we have that M* is a polyhedron with k facets, which 
we denote l-cosimplex. 

Further for M d S"^ convex we have the relations 



In this terminology we have for an instance A of F' and set Q = {ctji, . . . , 
ttj, } from step 5 of Algorithm |2.6| with Pa ^ ^ that M = co Q is an / — 1- 
simplex with M G PX and cos(i?'(M)) = def Q. Correspondingly M* is a 
/ - l-cosimplex with Pa C M* and t'IPa) < r'{M*) = % - R'{M). 

Moreover for M G S'^ the quantities R'{M) and D{M) are not indepen- 
dent as the following theorem by L. A. Santalo [jll| shows: 

Theorem 3.1 (Santalo) Let M g S'^ he contained in a half-sphere. Then 
1. For cos R' > 1/VdTl 



D{M) = sup{Z{x,y) \x,y e M}. 



M* = {xe S'^ \ m^x > for all m G M}. 



r'(M) + R'{M*) 



d{M) + D{M*) = 2n. 



cos2i?' < cosD < 



{d+ 1) cos' 



,^R'-1 



d 



2. For < cosi?' < l/^/d+l 



cos 2/2' < cos-D < 



{d + 1) cos^ R' - 1 
1 + (rf + 1) cos2 R' 



for d odd, 



cos2i?' < cosD < 



{d + 1) cos' 



•?R' -I 




(rf+l)(d-2) 
d+2 



cos2 R'J 
for d even. 
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Using Santalo's theorem we may estimate the breadth d{PA) of Pa by 
computing R'{Q). After a suitable transformation of coordinates we have for 
all (xi, . . . ,Xd+i)'^ € Pa that l^il < sin d{PA), i.e. all points of Pa are close 
to the equator {(xi, . . . , Xd+i) E S'^ \ Xi = 0}, and Pa can be enlarged by a 
transformation. If we know that in the beginning Pa is not to small, we must 
find a feasible x'' after a number of steps, for otherwise the transformed Pa 
would not fit on the sphere. 

A suitable quantity to measure the size of Pa is the volume and the es- 
timate of the size of Pa at the start of the algorithm is standard from the 
theory of linear programming. This leaves the determination of the hyper- 
plane. While the proof of Santalo's theorem is constructive, the construc- 
tion of the hyperplane needs the computation of about (^^2) distances and 
is therefore not practicable for a polynomial algorithm. Thus we substitute 
Santalo's theorem by a weaker one, for which we can easily compute all 
quantities. 

Definition 3.2 Let oi, . . . , a^+i G 5*^, such that Q = co {oi, . . . , a^+i} is a 
spherical simplex. For each at let a[ he the intersection of the circle through 
Qi and the circumcenter C {Q) with co {ai, . . . , aj_i, aj+i, . . . , a^+i}. Then the 
vertex diameter Dy{Q) is defined by 

D^{Q) = max{Z(ai, a[),..., /(a^+i, a^+i)}. 

Clearly the vertex diameter can be computed easily. While contains 
no explicit estimate of the vertex diameter, the methods can be adopted. To 
be complete we present the changes in the proof. First we have for the vertex 
diameter of the regular simplex 

Lemma 3.3 (Santalo) Let Q G S'^ be a regular simplex with circumradius 
R' and vertex diameter D^, then 

jd + l) cos^ R' - 1 

(l + (rf-l)(rf+l)cOs2 ^ ' 

Proof. The formula is a simple combination of formulas (2.3) and (2.13) in 
ilTI. □ 



Lemma 3.4 (Santalo) Let Q = co {ai, . . . , 0^+1} C S"^ be a spherical sim- 
plex, such that all vertices are contained in its circumsphere. Let C = 
Yl'i=i f'i'^i ^6 the circumcenter of Q. Then for the circumradius R' of Q 
it holds 

(d+l \ "1 ^ d+l 

i=l / i=l 
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Proof. We have 

d+l d+l 



cosR' = ^/iiofC = C^C = 1. 



i=l i=l 

From the equation we obtain 

d+l 



i=l l<«<j<d+l 

The inequahty is an immediate consequence of 2/ij/ij < + ij?,. □ 



Lemma 3.5 (Santalo) Among all spherical d-simplices on S"^ with circum- 
radius R' and cos R' < 1/ \/d+ I, which contain all vertices in their circum- 
sphere, the regular simplex has minimal vertex diameter. 

Proof. Assume that there is a simplex Q G S"^ with circumradius R', circum- 
center C and vertex diameter A„, which has all vertices on its circumsphere 
and has smaller vertex diameter than the regular simplex with circumradius 
R'. Let be the vertex diameter of the regular simplex with circumradius 
R'. From (0) we have cos-D^, < 0. Let ai, . . . , aa+i be the vertices of Q. We 
have 

d+l 

C = ^^yUjOj with /ij > 0. 

2=1 



and 

From this we obtain 
Hi cos R' = fiiaJC 



d+l d+l 
i=2 1=2 



'd+l N 1/2 



= /ii + /ii I ^ yU- + 2 ^ fiifijafaj j cos Z(ai, a[). 

\i=2 2<i<j<d+l J 

By our assumption we have cosZ(ai,a'^) > cosD^. Thus yields 

/d+l \ 

/ii cos i?' — /i^ — /ii j /ij + 2 fiifijttjaj] cosD^ > 0. (9) 

\i=2 2<i<jr<d+l / 
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For 02, ... , aa+i we obtain inequalities analogous to (0). Summing up all 



these inequalities, using the Cauchy-Schwarz inequality x y < \\x 
and taking ||C|| = 1 into account we obtain 



|l/2|U,||l/2 



d+1 



d+1 



d+1 



( 



1/2 



< ^/ifccosi?' -^^l -^^k 



k=l 

d+1 



k=l 
d+1 



d+1 



T 



i=l 
\i^k 



l<i<j<d+l 
ij'/'k 



COS Dy 



< Hk COS R' — 

k=l 

/d+1 \ 1/2 



l4 



k=l 
d+1 



i+1 



1/2 



cos A) 



\A;=1 

d+1 



i=l 
d+1 



i=l 



'd+1 



1/2 



^ fik COS R' -^fil- 



k=l 



k=l 



.k=l 



/d+1 \ 

I^IJ''! + (d-l)] cosDy 



ma 



(10) 

Solving (|l^) for cosD^ and using the identity and inequality from Lem- 



3^ gives finally 

cos A, < 



< 



{d + 1) cos2 R' - 1 



1/2 



(1 + l)(rf- l)cos2 R'Y^^ 
This is a contradiction to (^. 



□ 



By Lemmas |3]^, |3.5| we may determine with the help of a set Q two half- 
spheres, which contain P4 in their intersection and the angle (p between them 
is most 

l-{d+l) cos2 R'{Q) 



cos 6 > 



(11) 



l + (d-l)(rf+l)cOs2 /?'(g))l/2' 

(pUD yields by an easy calculation and further estimation 

sin0 < (d + l)cosi?'. (12) 

Next we show that for sufficiently small we may find a transformation 
of the sphere, which maps spherical polyhedra into spherical polyhedra and 
increases the volume of Pa by a prescribed factor. Any spherical map 

B-.S'^^S'^, B{x) = B'x/\\B'x\\ 
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with a nondegenerate linear map B' has the first property. The influence of 
B on the volume is most easily studied by use of parametrisations. 

A parameterized rf-surface in R'^+i is a difTerentiable injective map $ : 
G C R"' ^ R°'+-'^ such that for every x E G the Jacobi-matrix {^) is 
nondegenerate. For M C ^{G) the volume of M is given by 

V{M) = / det $ dx, 

where det $ denotes the volume of the d-dimensional parallelepiped spanned 
by . . . , For a differentiable map ip : R'^+-^ — > with appropriate 
properties ^ o $ is a parameterized surface and consequently 



V{^{M)) = I det(^o$)da; 



'*-i(M) 

is the volume of ^^(M). Thus (det(^' o $)/ det $)($"^(a;)) is the local change 
of volume in x e $(G) under the map ^ : ^ ^''"'"^ 

Here we are interested in the case 



: 1R'=^+^ - 



l-Y^xif, (13) 



(14) 



'^a{xi, Xd+i) = {axi, X2,..., Xd+if / \Jl + {a^ - l)xl, a > 0. 



The restriction of on S'^ is the spherical transformation {axi,X2-, ■ ■ ■ , 
^d+i)'^ /\\{'^Xi,X2, ■ ■ ■ ,XiiJ^_i)'^\\. We have chosen rather than the other 
transformation, as its derivatives can be computed more easily. In particular 
we consider the point x — (/?, 0, . . . , 0, \/l ~ I3'^Y ^i^^ < /3 < 1. Compu- 
tation of the Jacobi matrices and of the determinants with the help of Gram 
determinants yields 

1 a 

det $(x) = — — , det(*e«o$)(^) = 



(1-/32)1/2' ^ " ' (l-/32)l/2(l + /32(Q,2_l))(<i+l)/2- 

Thus the transformation tl^a gives for all {x\, . . . ,Xd+i f e S'^ with xi = 13 
the local change of volume 

m 



l + /?2(a2_l))(d+l)/2' 
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Consequently we obtain for M (Z S"^ with x\ < (3 for all (xi, . . . , Xd+i Y EM 
and a > the estimate 

We may summarize the previous considerations in 

Lemma 3.6 Let M G S'^ such that for all {xi, . . . ,Xd+i)'^ E M it holds 
Xi < (3. Further let '■ S'^ S'^, a > 0, be given by (|7^. Then the volume 
o/\I'q,(M) can be estimated by j^ldj). 

Evidently we may generalize the transformation to the case of a trans- 
formation orthogonal to an arbitrary equator a'^x = instead of the equator 
xi = 0. 

As we work with the homogenized problem, we have to study the influence 
of on the volume. Using the same technique as before we find for M C M*^ 
for the volume of the homogenized set 

Thus we have 

Lemma 3.7 Let M C R'^ with \\x\\ < R for all x E M. Then V{^{M)) > 

The following considerations concerning polynomial algorithms could be 
done in the context of oracle-polynomial time algorithms for convex bodies in 
The presentation becomes somewhat simpler if we follow the presentation 
in for algorithms for the feasibility problem of linear systems. 

According to the corollary to Lemma 8.7 in P| we can solve F in polyno- 
mial time, if we can solve the following problem in polynomial time: For the 
set 

M = {x eR'^ \ ajx > hi, i = 1,. . .,n} (16) 

of size L find a point in M or show that there exists no such point in time 
polynomial in L. Here the size L of the system measures the amount of 
space needed to write down the inequalities. For a precise definition of the 
size compare H| or [^. 

We have that either M is empty or the volume of M is not too small: 

Lemma 3.8 // the set M from ( ff^ j is nonempty, then 
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Proof. 0, Lemma 8.14. □ 

For the a^, bi from (|16D let a[ = {aj , —bi)'^. Then we have for the homog- 
enization of M 

$(M) = {x e S'^l afx > 0, i = l,...,n, Xd+i > 0}. 
We may combine Lemmas |3.71 , |3.8| to obtain 
Lemma 3.9 Let M be the set from I ^Td{ ) with size L. Then 

l^($(M)) > 2-3('^+2)L 

or $(M) is empty. 

Now we have collected the parts to show that the following modification 
of Algorithm |2]^ is polynomial. 



Algorithm 3.10 Let (3^ = V((4/3)2/('^+2) - l)/3, max = 6{d + 2)L. 
Main loop: 

1. Letp = 0. 

2. Let p = p + 1. If p > max, stop with infeasibility. Let j G {1, . . . , n} 



be arbitrary, k = 1, x^ = aj , Qi 



3. If x^/Wx'^W is feasible, then stop with feasibility . Else let m be the index 
of any constraint, which is violated for a;*^/||a;'^|| . // Qk U {am} is a 
positively spanning set, then stop with infeasibility. Let y be the point 
on the line segment x^am, which is closest to the origin. 

4. Compute the center C of the touching sphere of Qk U {cim}- 

5. If C E conv {Qk U {am}), then do the following: Let x'^^^ = C , Qk+i = 
Qk U {am}, k = k + 1. If def (Qk+i) < Pd/{d + 1) then Transform and 
goto step 2 else goto step 3. 

6. Let y be the intersection of the line segment yC with the relative bound- 
ary of cony {Qk U {am})- Let F be a facet of conv {Qk U {am}), which 
contains y. Let aj be the vertex of conv {Qk U {am}), which is not 
contained in F. Let Qk = Qk \ {c-j}- Goto step 4- 

Transform: 

1. Determine the vertex diameter of co Qk. Let ai be the vertex which 
gives the vertex diameter. 
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2. Apply the transformation \&2 with respect to the equator afx = 0. 



Theorem 3.11 Algorithm solves F in polynomial time. 



Proof. We denote by Mp the set of feasible points before the p-th transfor- 
mation. First we show the correctness of the algorithm. The first stopping 
condition in step 3 gives trivially the correct result. If the algorithm stops 
at the second stopping condition, then M C Pq^,. Pq^. is a lower dimensional 
sphere and thus ^(Pqj.) = and the instance is infeasible by Lemma pIOl . 

Now we consider the state of the algorithm at the call of the procedure 
"Transform". We have cos R'{Qk) < f3d/{d+ 1). By (pD we find that all 
points of Mp are contained in the zone {x E S'^ \ < ajx < Pd} for the 
Oj giving the vertex diameter. From this we obtain by Lemma |3]^ that 
V{Mp^i) /V{Mp) > 3/2. Finally we obtain by Lemma ^]9| that after max 
steps we have V{Mmax) > 2 > V{S'^). Consequently an existing feasible 
point must be found before this step. This proves the correctness of the 
algorithm. 

To calculate the number of computations we only have to estimate the 
number of iterations in the inner loop starting in step 2. By a power series 
expansion we see easily that P > cd'^^"^ for suitable c > 0. Thus by Lemma 
|2.11| there are 0{d^) iterations. 

It remains to show that it is sufficient to do all calculations with rational 
numbers whose size is bounded by a polynomial in L. This can be done in 
the same way as in Khachiyan's algorithm and we skip the details. □ 

Some remarks concerning the algorithm: 



1. If we compare Algorithm p.9| and Algorithm |3.10| we find that we have 
replaced the condition in Algorithm |2]^ that the new in step 2 is 
given by a most violated constraint by the weaker condition that the 



new Qm in step 4 of Algorithm |3.1CI| is given by any violated constraint. 



We can easily check that it is sufficient to choose any a 7^ 0, such that 
a^x^ < 0, but a^x > for all x G Pa, i.e. if {a^x = 0} is a separating 
great sphere for x^ and Pa- Moreover we made no use of the fact that 
Pa is a polyhedron, but used only that the volume of Pa cannot be 
too small. Now a separating hyperplane for the unhomogenized prob- 
lem immediately yields a separating great sphere for the homogenized 
problem. Thus the algorithm is essentially a polynomial algorithm for 
the following convex feasibility problem: Given is a convex set K & 
which is either empty or well bounded in the language of 0] and an 
oracle which for any a; G M*^ confirms that x & K 01 gives a separating 
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hj^erplane. Under these assumptions find a point in K or show that 
K is empty. 



2. While the previous remark pointed out that from a theoretical point 
of view, it is sufficient to choose any separating great sphere in step 
4, from a practical point of view this is certainly not the case. This 
is shown by the following example, which shows that the inner loop in 
Algorithm |3.1CI| may need many iterations: 

In this example the sets Qk have always 3 elements which are for sim- 
plicity unnormalized. Let 

Qfc = {(1, e, 0, . . . , 0)^, (-1, e, 0, . . . , 0)^, (0, -4, 1, 0, ... , 0)^}, 

where e and 6k are positive numbers of the kind that 6k is small with 
respect to 1, and e is small with respect to 6k. The points are of 
the form x'' = (0, 1, r^fc, 0, . . . , 0)^. {(0, -6k+i, 1,0,..., 0)^x = 0} is a 
separating great sphere if 6k+i > rjk- The computation of x'^^'^ gives 



+ 

= (0, 1, Tik+i, 0, . . . , 0)"^, r/fc+i = 4+1 + e \ . (17) 

Vl + 

Thus if 4+1 is close to r]k we have that the distance /^{x^,x^^^) is of 
order e. Now the iteration may be repeated with a 4+2 which is of the 
same order as 4+i- This behavior does not affect the polynomiality 
as e is bounded below by a function depending on the dimension, but 
demonstrates that the convergence may be rather slow. 

3. If we choose in step 4 the most violated constraint then this phe- 
nomenon cannot occur, as is shown by the discussion after Lemma 



2.11 



4. While the choice of x^ as circumcenter of the sets Qk is fundamental 
for our work, this is certainly not the case for the construction of the 
transformation. Although it did secure the polynomiality of the algo- 
rithm, there are some apparent disadvantages. We have to work with 
three different quantities (circumradius, vertex diameter, volume) and 
these quantities are rather weakly related. Moreover the volume does 
not behave very nicely for transformations of the sphere. The apparent 
choice for a transformation would be to transform after each iteration 
coQk into a regular simplex of the same circumradius such that the cir- 
cumcenter of CO Qk is transformed to the circumcenter of the resulting 
regular simplex. Unfortunately we don't know of any quantity which 
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behaves nicely under this transformation and thus don't even know, 
whether the resuhing algorithm is convergent. 



5. Still there are probably other transformations which are easier to handle 
than the one sketched above and lead to polynomial algorithms. One 
approach might be the following: We want to construct in polynomial 
time a positively spanning set for infeasible problems. For our heuristic 
considerations we assume that this positively spanning set is actually a 
positive basis. As we want to construct a positive basis, we look at the 
violation. Considering the fact that for a regular simplex the quotient 
of the circumradius and the inradius is d we find that f (x'^) > l/((i+l), 
if the points of the positive basis form a regular simplex. Now, let us 
assume that we have found a point such that v{x^) is small compared 
to + 1). Let Q be an arbitrary positive basis spanned by the a,. 
Then it is easy to see, that there must be a face F of convQ, such 
that F is nearly parallel to {{x^Y'x = 0} and has a small distance from 
the origin. Let us apply for suitable a > 1 a transformation with 
respect to the equator {{x'')'^x = 0}. Locally we increase the violation. 
Globally it may be seen, that the distance to the origin of all faces of 
convQ, which have a dimension less or equal than the dimension of F 
is increased. One may hope that the minimal rate of increase is given 
by a function, which depends only on the dimension. Thus one should 
obtain fast, i.e. polynomial convergence. 

We have not tried to work out the theoretical details of this algorithm. 
Thus we don't know whether it works and have no information on 
its theoretical quality. We may say that this algorithm would avoid 
the problems of our algorithm which we pointed out in the previous 
remark, as in our context it is certainly more natural to deal with dis- 
tances than with volume. Moreover the transformation with a pole 
is more natural than the transformation with a pole Oj and it is easy to 
implement this algorithm. This was in fact done with a transformation 
taking place whenever v{x'') < 1/^/d. The surprisingly good results of 
this algorithm are presented in the next section. 

Altogether we may say, that solely as a consequence of the weak estimate 



in Lemma p.llj it is possible to construct a sequence of transformations which 
transform the combinatorial Algorithm |2.9| into a polynomial one. This is a 
remarkable contrast to the combinatorial simplex algorithm, as for the sim- 
plex algorithm it is even unknown, whether an optimal choice of pivot-steps 
would lead to a polynomial algorithm (HiRSCH-conjecture, for an account 
see PI). 
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Moreover there is some hope, that Algorithm ^]6| is polynomial. That 
would be of particular interest, as an algorithm which is combinatorial as 
well as polynomial might be a candidate for a strongly polynomial algorithm 
for F. It is still an open question whether such an algorithm exists. For 
details about this problem see 0. 



4 Implementation of the algorithms and nu- 
merical experiments 

Up to now we have only studied the geometry of the algorithms and we still 
have the problem to perform the necessary calculations. We don't want to 
go into details or find the optimal way to do this but just show that there are 
no particular problems and that the complexity of a single step is on average 
comparable to the simplex algorithm, i.e. 0{dn). 

The only steps in Algorithm ^]6| where the computations are not com- 
pletely obvious are step 3 and step 4. Here we must compute the center C of 
the touching sphere and check whether it is in the convex hull of the Oj. We 
observe that by our construction is linearly independent and Qk U {am} 
affinely independent. We compute a least squares representation 

Qm = a + Xitti, a orthogonal to all G Qk (18) 

and distinguish two cases. 

1. a = 0. Then Qk U {am} is linearly dependent, C = and using (p^) we 
easily find a representation of as an affine combination of Qk U {am}- 

2. a 7^ 0. We observe that C is characterized by the following three 
properties: 

(a) aJC = aJC > for all a^, aj G Qfc U {am}. 

(b) C = Ea,eQ,u{a„} with m e M. 

(c) E/^i = 1- 

Thus if a C" 7^ has the first two properties, then C is just a multiple 
of C. 
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Let us look at the points y constructed in step 6 of the algorithm. We 
find inductively that we have 



Z(?/, ai) = Z{y, ai) < Z{y, a„) for all a^, ai G Qk 
y = fXitti, with all /ij > 0. 

aiGQfcU{am} 

Thus using ( ]T8| ) we find that there is a C of the form C = y + pa for 
suitable p > 0. From this we determine C as an affine combination of 

Qk- 

The next aj to leave Qk is then apparently found by a simple line search. 
It remains to solve ([T8|). To do this we consider the Qk as matrices with 
columns Oj. We observe that consecutive matrices Qk, Qk+i only differ in 
that a column is added or taken away. If we solve (0) by computing a 
QR-factorization of Qk then it is well known that the QR-factorizations of 
consecutive Qk, Qk+i can be obtained by an updating process which takes 
only 0{d?) computations, cf. [Q. Thus taking into account that the deter- 
mination of in step 2 takes 0{dn) computations and we have on average 
to add one column to Qk and take one column away, we find that a step of 
Algorithm |2.6| takes 0{dn) computations on average. 

Clearly our results above give no information whether the algorithms 
could be of practical use. An obvious way to test this would be to solve the 
test problems as provided by |http: / / www.netlib.org/lp| . These problems are 



sparse optimization problems, i.e. to solve them we have to deal with many 
issues which are rather unrelated to the basic algorithm, but will strongly in- 
fluence the results. Among these are solving the optimization problem by an 
algorithm for feasibility, dealing with equality-constraints, dealing with the 
sparsity. Thus we have taken a much more simple approach. We generated 
three different sets of problems. The most basic way was done by generat- 
ing in a random way n vectors of unit length and n negative constants 
bi (Ex 1 in the following tables). Thus we obtain a system of inequalities 
afx > bi which has the origin as a feasible point. In the next problem we 
chose 6j = for i = 1, . . . , (i + 1 and a^+i = — Yl'i=i '^i/W Yl'i=i Thus in 
this case we obtained with probability 1 problems which have exactly one 
feasible point (Ex 2). The last set differed from the second only in that we 
chose bd+i = 10 to construct infeasible problems (Ex 3). In order to use the 
origin as a starting point all systems were translated by a random vector. 

If we look at the problems then we see that after choosing the correct 
set of indices the second problem is equivalent to the solution of a system of 
linear equations. If we use the Gauss-algorithm, and allow to transform the 



20 



complete problem in every step then we see that a solution would involve d 
steps, each with 0{nd) operations. In particular the number of steps would 
be independent of the number of inequalities. Certainly we would consider 
such an algorithm as excellent from a practical point of view and we may ask 
how close we can get. 

We used three different algorithms for the solution of the problems. First 
we used the following variant of the simplex algorithm to have a basis of 
comparison. We consider the problem as the dual of a problem in standard 
form. We perform d pivot-steps to generate a basis. This corresponds to ask 
for equality in the inequalities corresponding to the basic columns. Then we 
choose a positive constant right side and use the simplex-algorithm to solve 
the generated problem. If the algorithm stops with a finite solution, this 
shows the existence of a feasible point in the dual problem in the intersection 
of the constraints corresponding to the basic columns. If the algorithm stops 
with unboundedness this shows the infeasibility of the dual problem. This 
algorithm is denoted Alg 1 in the following in the tables. 

The next algorithm was a straightforward implementation of Algorithm 
p.6| (Alg 2). The last experimental algorithm (Alg 3) used an additional 
rescaling in Algorithm p.6| . Here it does not make much sense to use the 
polynomial Algorithm |3.1CI| as the number of steps before the first transfor- 
mation in step 5 takes place is not much smaller than the total number of 



steps for the complete solution by Algorithm 2.6. Instead we used the fol- 



lowing procedure. We have by (|^) that every step of the algorithm makes 
relative progress of at least ■>/! — f (x^)^. Thus when f (x^) becomes less than 
l/\/d then we do the following. Let v{x'') be determined by the constraint 
Oj., i.e. f(a;'^) = — a^x'^/||x'^||. We set a- = {I + Xx^x^'^)ai, i = l,...,n 
and renormalize the a'- to unit length. Here A is chosen in a way, such that 
a^x'^/llx'^ll = —y^2/d. We observe that the matrix Qk is transformed by 
the formula Qk = /i(<5fc + wx^'^) for a scalar /i and a vector w which are 
easily determined. For such a transformation it is possible to update the 
QR-factorization in 0{d^) computations, cf |Q. Thus the transformation 
takes a total of 0{nd) computations. 

To test the dependence on the dimension we chose the dimensions d = 10, 
20, 40, 80, 160, 320 and 640 and 8d inequalities for each dimension. To 
test the dependence on the number of inequalities we chose d = 100 and 
n = 400, 800, 1600, 3200, 6400. For each case 5 examples were generated and 
solved. The average numbers of steps are shown in the tables. 

In the tables, the second column for the third algorithm gives the number 
of rescalings. Moreover the numbers for the simplex algorithm do not give the 
d steps necessary to obtain the starting tableau. We see that all algorithms 
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Table 1: Number of steps in dependence of the dimension 

show httle dependence on the number of constraints, in particular for Alg 3 
this dependence appears to be sub-logarithmic. 

For the dependence on d we may not expect spectacular results, as we 
know by the remarks at the beginning, that we may not hope for less than d 
iterations and on the other hand the simplex algorithm solved all examples 
in less than 12d iterations. But we see that for large dimensions the new 
algorithms need fewer iterations and in particular Alg 3, the algorithm with 
the rescaling, did very well. 

To obtain some more insight in the dependence on the dimension we did 
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1.1122 1.3093 
1.4575 1.2719 
1.3538 1.2747 


1.6269 1.0214 
2.4228 1.0334 
2.2439 1.0334 



Table 2: Rate of growth 
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Table 3: Number of steps in dependence of the number of inequalities 



a best fit of the form ad'^ for all series and all algorithms. Here we computed 
two values for the simplex algorithm, for the first value we counted the d steps 
for the initialization while in the second value these steps were neglected. For 
the algorithm with the rescaling we did not take into account the numbers of 
rescaling, as this number appears to have little dependence on the dimension. 
The results are presented in the second table. 

We find that the simplex algorithm and Algorithm p.6| are superlinear, 
where the exponent of Algorithm ^]6| is slightly smaller for all series, while 
the experimental algorithm is very close to being linear for our examples and 
thus for these examples very close to the hypothetical Gauss-type algorithm. 

Altogether the examples show that the concepts developed here might be 
of practical interest. 

Finally the author thanks M. Henk and A. Schurmann for valuable 
discussions, hints to the literature and reading a previous version of this 
paper. 
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